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A simple numerical model is used to simulate the effect of vertical taps on a packing of monodis- 
perse hard spheres. Our results are in good agreement with an experimental work done in Chicago 
and with other previous models, especially concerning the dynamics of the compaction, the influence 
of the excitation strength on the compaction efficiency, and some ageing effects. The principal asset 
of the model is that it allows a local analysis of the packings. Vertical and transverse density profiles 
are used as well as size and volume distributions of the pores. An interesting result concerns the ap- 
pearance of a vertical gradient in the density profiles during compaction. Furthermore, the volume 
distribution of the pores suggests that the smallest pores, ranging in size between a tetrahedral and 
an octahedral site, are not strongly affected by the tapping process, in contrast to the largest pores 
which are more sensitive to the compaction of the packing. 

PACS numbers: 45.05.x, 45.70.Cc, 61.20.-p, 81.05.Rm, 81.40. Cd 



I. INTRODUCTION 

Granular materials constitute the raw materials in a 
huge number of human activities as in agriculture, the 
mining industry, pharmaceutics and are at the heart of 
the matter in several ecological concerns as desertification 
by eolian erosion or avalanches. Therefore, explaining a 
few current granular processes, such as storage, trans- 
port, or collapse, is a real economical challenge. Further- 
more, packings of spheres which is the simplest model 
for granular medium, have a great fundamental interest 
for physicists: hard sphere systems are indeed a common 
description of simple liquids jlj ; moreover grains can be- 
have, according to the external conditions, more or less 
like a solid, a liquid or a gas This great variety of 
behaviors for a banal heap of grains makes granular me- 
chanics a rich area of investigation only partially clarified 
at the moment. It is now a well-known result ||-^| (al- 
though there is no theoretical explanation for it) that a 
disordered static packing of equal hard spheres can cover 
a large range of volume fraction, approximately from 56 
%, the random loose packing (R.L.P.), to 64 %, the ran- 
dom close packing (R.C.P.). For a regular arrangement, 
the packing fraction can reach up to 74 % which cor- 
responds to the densest structures, namely the hexago- 
nal compact (H.C.) and the face-centered cubic (F.C.C.) 
crystals. 

As thermal energy (fcgT) plays no role, because it 
is insignificant compared to the gravitational energy 
of a macroscopic grain, each packing of spheres is a 
metastable configuration which can persist as long as 
there is no external excitation. In this frame, the is- 
sues of compaction of grains under vertical taps are a 
practical way to study the succession of jumps from a 
metastable equilibrium to another one. The initial pack- 
ing is quite loose and can progressively reach a nearly 
stationary configuration (steady state) evaluated through 
its average volume fraction. Some experiments done in 



Chicago [|]-|| have studied the influence of the tapping 
intensity on the steady-state value and the dynamics of 
the compaction, which is approximately inverse of the 
logarithm of the number of taps. The experimental set- 
up is a thin tube of diameter D=1.88 cm filled to about 
an 80 cm height with monodisperse, spherical soda-lime 
glass beads (of diameter d=l, 2 or 3 mm). The tube is 
shaken by an electromagnetic exciter delivering vertical 
taps, each of them consisting of an entire cycle of a sine- 
wave of frequency f = 30 Hz. The excitation strength 
is parameterized by T, the ratio between the measured 
acceleration peak and the gravitational acceleration g. 
Moreover, several numerical and theoretical works |)| 15] , 
most of them dealing with the notions of free volume 
and geometric constraint, found the same kind of behav- 
ior as obtained experimentally and some of them pO 11 
have pointed out structural ageing effects, as typically 
observed in glassy systems. So a parallel might exist 
between this granular compaction and the dynamics of 
out-of-equilibrium systems as glasses. 
In this work, we used a simple model to simulate the 
compaction of a packing of monosize spheres submit- 
ted to vertical taps. We did not try to make a realistic 
description of the quite complex succession of collisions 
in a shaken packing: we have kept as the only ingredi- 
ent of the model the geometric constraint between hard 
spheres, which is believed to be the principal origin of the 
compaction. Despite the fact that we deliberately forgot 
the mechanical dimension of the problem, the model is 
able to reproduce qualitatively the experimental results 
of the Chicago group and some further results in agree- 
ment with different numerical and theoretical studies. As 
the model seems to capture the physics of the problem, it 
is then possible to go beyond a global analysis. Indeed, as 
a three-dimensional packing of hard spheres, our descrip- 
tion has the quite interesting asset that it is very close to 
a real granular medium. So, contrary to almost all the 
previous works which deal only with a macroscopic probe 
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(i.e. the average density in all or part of the packing), 
our model can provide us with realistic information on 
the local structure of a packing and its evolution under 
compaction by taps. 

The paper is organized as follows. A detailed description 
of the model is presented in section II. Section III is de- 
voted to the global analysis of compaction (logarithmic 
dynamics, hysteresis effect, and ageing behaviors). In 
section IV, the local analysis of the packings is described 
with the use of density profiles and size and volume dis- 
tributions of the pores. Our conclusions and perspectives 
end the paper in section V. 



II. THE MODEL 

The model proposed here is purely geometric and deals 
only with the steric constraint, neither friction nor con- 
tact law between the spheres or with the walls is intro- 
duced. 

The different sequences of tapping were initialized from a 
relatively loose packing obtained by a steepest-descent al- 
gorithm simulating a sequential gravitational deposition 
[M . We worked with packings of 4096 spheres of radius 
R piled up in a square-box of dimension L = 32R. Con- 
cerning the vertical walls, we used both periodic bound- 
ary conditions (P. B.C.) and fixed boundary conditions 
(F.B.C.) i.e. impassable vertical planes. The top of the 
box is open whereas the bottom is a fixed impassable 
plane. 

A tap is decomposed in two stages: first a vertical dila- 
tion and then a gravitational redeposition. 
The first stage corresponds to the external excitation 
which will enable the packing to move from a metastable 
equilibrium to another one. We used the simplest way to 
simulate the tap by applying an uniform dilation e to the 
whole packing (z — > z(l+e)). This reduction is certainly 
far from a real tap but we assume that the way of dilat- 
ing the packing is less important than the result of the 
dilation: a significant increase of the average free- volume 
of the spheres which will allow collective rearrangements 
during the second stage, the redeposition of the packing. 
This redeposition procedure must be non-sequential in 
order to permit such collective behaviors; so we use a 
monte-carlo algorithm to discretise the motion of the 
spheres: a great number of small displacements are com- 
puted . An individual movement procedure is structured 
as follows: a sphere, randomly chosen, is submitted to 
a small random displacement; if this displacement cre- 
ates no interpenetration with another sphere or with the 
walls (according to the boundary conditions), it is ac- 
cepted otherwise it is rejected. Because of this binary 
schema, two neighboring spheres can not be exactly in 
contact but, after a sufficient time, they get very close 
to contact. Figure |l| shows a typical displacement: the 
values of the polar angle <fi and of the displacement d 
are strictly randomly chosen respectively between and 



2ir and between and d max , whereas the choice of the 
angle 9 follows a random distribution centered on zero 
to mimic the effect of gravity. We used the following 
Gaussian distribution of width 0q truncated beyond 7r/2 
in order to orientate all the displacements down to the 
bottom of the box. 

P(6)=AcM-(0/0o) 2 ) (1) 

The choice of the distribution docs not seem to be re- 
strictive: some attempts with a Poissonian and a linear 
distribution give qualitatively the same phenomenology; 
the pertinent parameter is the width 9q. 
With such an algorithm, agitation will persist indef- 
initely So we regularly test the packing during the 
redeposition process. The variable checked is (Z), the 
average altitude of the packing that is the average po- 
tential energy of the spheres. The redeposition is stopped 
when the relative variation of (Z) becomes smaller than 
a threshold 77. The choice of (Z) is motivated by its easy 
evaluation during the process and by its possible link 
with a statistical mechanics approach. 
This simulation is rather close to the one proposed by 
Barker and Mehta but with some differences es- 

pecially concerning the way of introducing gravity and 
the end of the redeposition stage. 

The model uses four parameters: d max , 77, 9q and e. 
The two first ones have a direct effect on the simulation 
time. The smaller 77, the longer the simulation time; 
still, 77 must be small enough if we want the redeposition 
to be nearly completed. The parameter d max has to be 
optimised. A very small value of d max allows almost all 
of the displacements to be accepted but the effect on the 
redeposition is very slight and the packing is therefore 
nearly frozen. On the contrary, for a large value of d max , 
almost of all the displacements are refused and, once 
again, the packing evolves very slowly. In this study, we 
used the intermediate value d max = R/5. 
80 has a significant effect on the packing behavior: a very 
small Oo induces a decompaction whereas a large value 
decreases the efficiency of the compaction. We found 80 
— 7r/4 as the optimised value giving rise to the maximal 
compaction rate. 

The last parameter, e, corresponds to the external ex- 
citation induced in the packing. This is our control 
parameter. The value of e can be estimated from ex- 
perimental results concerning the dilation of a vertically 
shaken sand heap @: £ = Jh/h « 5/500 « I0~ 2 . We 
can also try to link roughly e to the experimental control 
parameter, the dimensionless acceleration T — Aw 2 /g 
where A and ui are, respectively, the amplitude and the 
frequency imposed to the bottom of the heap. In first 
approximation, if we neglect the loss of energy in the 
packing, a particle at the top of the heap (z(0) = h) ac- 
quires an initial speed wA and achieves a ballistic flight. 
Its maximal altitude is z(0 + ) = h + | (^) 2 and then it 
comes: 
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T = u;( — ) 1 /V/ 2 =► Tcxe 1 / 2 
9 



(2) 



X(t) = X oc - 



As e 1 / 2 is linked to T, we will use it as our control pa- 
rameter to quantify the strength of the tapping process. 
With this, it is possible to compare the results of our 
model with the experimental work of the Chicago group 
and with other numerical and theoretical models, almost 
all of them dealing only with a global description of the 
granular system. 



III. GLOBAL ANALYSIS 

This global analysis is achieved with different average 
values. We did not use a direct evaluation of the packing 
fraction from the number of spheres in a reference volume 
because wether boundary effects are significant or, for a 
smaller volume, the statistics become too poor. More- 
over, the choice of the reference volume is not unique: 
it can be, for example, the space that contains all the 
spheres or the smaller one that contains only the cen- 
ters of the spheres. To avoid being partial, we evalu- 
ate the packing fraction by averaging the surface packing 
fraction, <I>, calculated on many horizontal cuts. This 
measure is permissible because of the following stereo- 
logic result: the average surface fraction of any cut in a 
pac king is equal to the volume fraction of the packing 
p0[ ; with horizontal cuts, this calculus is just a spatial 
integration which gives the exact volume fraction. The 
quantity $b is calculated in this way at the bottom of the 
packing between the heights and 4R; ($ c ) comes from a 
similar calculation on approximately 90 % of the packing 
and is corrected near the bottom wall by a perturbated 
This model uses a corrective factor for 
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the average density of a packing near a wall (between 
and R) with regard of a packing not perturbated by any 
wall. For the case of spheres near a plane, this factor is 
estimated to 16/11. It is also interesting to study (Z), 
the average potential energy of the whole system, which 
is quite pertinent in a statistical mechanics description. 



A. The dynamics of compaction 

The densification of the packing is observed through 
the temporal evolution of the preceding mean values; 
here the time is the number of taps and what we call 
dynamics of the compaction is, in fact, the succession of 
metastable equilibrium, each jump from one to another 
being induced by the taps. Figure |^ shows compaction 
laws obtained with fixed boundaries (F.B.C.) and three 
different excitation rates. This excitation intensity e 1 / 2 
has a decide d effect on the compaction dynamics (see 
section III B ) . The simulation curves are in good agree- 



1 + B x ln(l + t/T X ) 



(3) 



with X = <!>(, or ( < I ,C ). For (Z), a nearly similar fit can be 
proposed: 



<Z)(t) = (Z) 



1 + B Z ln(l + t/r z ) 



<Z)oc 

<Z)o 



B z ln(l + t/r z ) 



(4) 



ment with the experimental data and compatible with 
the following fit previously proposed H : 



We have noticed that a sum of two exponentials can also 
fit (Z)(t) reasonably well. 

The dependance of these parameters on e is difficult to 
characterize. We simply note that the parameter B is 
consistent with an exponential dependance on e 1 / 2 (i.e. 

r). 

This compaction dynamics is quite particular: as the 
packing progressively densities, the compaction efficiency 
decreases. So the dynamics reduces speed and the system 
evolves to a steady-state without never really reaching 
it. This slowing down is particularly remarkable for the 
smallest values of e 1 / 2 . This specific dynamics requires 
the study of the densification on a logarithmic time scale. 
It is also interesting to analyse the fluctuations of the 
curves, especially when the packing becomes close to its 
asymptotic or steady-state (SS) limit. The power spec- 
trum of the fluctuations X-Xgs as a function of the fre- 
quency, i.e. the inverse of the taps number, shows more 
or less a simple power-law in a log-log diagram (with a 
slope in the range 1 to 1.5). The effect of e is noticable 
only for the high frequencies. Moreover, the simple stan- 
dard deviation of the fluctuations, ax = \/ ((X — X gg ) 2 ), 
seems to be directly proportional to e 1 / 2 or V. These re- 
sults, calculated for X = (Z), are presented in figure 0. 
Furthermore, we have noticed that the periodic boundary 
conditions do not qualitatively affect these observations; 
the same remark can be made concerning the results of 
the following section. 



B. Hysteresis on the steady-state values 

The next stage consists on studying the influence of the 
excitation parameter e 1 / 2 on the maximal value of the 
packing fraction. For this purpose, we carried out a suc- 
cession of simulations with a sequence of 4000 taps. The 
steady-state value is estimated by averaging the packing 
fraction on the 1000 last taps or directly through the last 
value. The smaller is e 1 / 2 and the larger is the difference 
between this steady-state value and the asymptotic value 
Xqq given by the fit. Moreover, in this small excitations 
range, Xoa can get over the R.C.P. limit and r, the char- 
acteristic time of the fit, increases spectacularly. In fact, 
the fit becomes more and more uncertain in so far as 
the triplet (X^, B, r) is no longer unique, and depends 
strongly on the range of taps over which the data fitting 
is performed. This deviation between the steady-state 
value and an uncertain asymptotic limit has also been 
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noticed in the experimental work of the Chicago group 
H and in some theoretical studies jlj|[l5|. The depen- 
dance of ($ c ) on e 1 / 2 is shown in figure |] (solid black 
squares). The different packings are obtained after 4000 
taps of strength e, starting from the same initial packing. 
The curve has a bell shape with a maximum between 0.1 
and 0.2. 

If we now compute a unique tapping sequence with a pro- 
gressive increase of the excitation e 1 / 2 after every 4000 
taps (constant excitation increment: Ae 1 / 2 = +0.025), 
we obtain nearly the same curve for (<1> C ) , as can be seen 
in figure ^| (open circles). When carrying out the same 
process in the opposite way i.e. with a progressive de- 
crease of e 1 / 2 (Ae 1//2 = —0.025), two things can happen: 
If, while increasing, e 1 / 2 went beyond a critical value of 
(e 1 / 2 )* w 0.15, the final packing fraction (<I> C ) does not 
decrease but increases a bit more to a maximum value. If 
we compute another increase process (Ae 1 / 2 = +0.025), 
we cover approximately the same values. This last upper 
branch, including the part above (e 1 / 2 )* is represented 
on figure | (up and down open triangles). As it is rela- 
tively well reproductible, it is called "reversible" . We can 
also notice on this reversible branch that ($ c ) decreases 
with e 1 / 2 . 

On the contrary, if e 1 / 2 stayed below (e 1 / 2 )* during 
the increase stage, the steady-state values do not evolve 
significantly; they are nearly frozen, and it is hard to 
estimate whether there is a compaction or a decom- 
paction process, because the dynamics is very slow. This 
last branch is called "irreversible" and reflects the great 
metastability of the corresponding packings. 
To summarize, there is a strong hysteresis effect which 
allows the maximum compaction rate to be reached by 
an e 1 / 2 increase-decrease sequence. These observations 
arc in very good agreement with the results of Nowak & 
al (t]]. In particular, Figure |] is to be compared to the 
experimental data obtained with 1mm diameter beads, 
corresponding to an aspect-ratio of nearly 19, close to 
that used in our simulation (L/2R = 16). Surprisingly, 
for an aspect-ratio of 9, the experimental results show 
a much larger increase of the packing fraction on the 
reversible branch, up to nearly 66 % (i.e. more than 
the R.C.P. limit which may indicate a commensurability 
between the cylinder and the beads jjj). However, for 
a still smaller aspect-ratio of 6, the reversible branch is 
more similar to the first case, with a moderate increase 
to a maximal value below the 64 % limit. 



and the same packing after an evolution time tw (waiting 
time) . In this study we work with the following function: 



A(t,^) = ((Z)(t)-(Z)(t + t w )) 2 



(5) 



Here, x indicates the statistical average of x; that is, the 
mean value calculated for a certain number of realizations 
of the same experiment. The results have been averaged 
on only 10 realizations because of the limitation due to 
the calculation time. The statistics are, therefore, rather 
poor, that is why we use solely (Z), which fluctuates quite 
less than the other global values. In figure are drawn 
the curves of A(t,tw) obtained for different values of tw- 
The re is o bviously a scaling law; a similar fit as in sec- 
. [II A , with the three parameters (the asymptotic 
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limit), Ba and ta, is quite compatible with the data: 



A(t,t w ) =A C 



1 + B A ln(l + t/T A ] 



(6) 



The same kind of ageing effects have already been pointed 
out in previous numerical studies p0| , pT| . These ef- 
fects confirms the great similarity between granular 
compaction or more generally slow granular rheology 
and glassy systems submitted to time-dependent driving 
forces (see for instance [f22|- 



To conclude with the global analysis of the compaction, 
it is satisfying to note that our simulation reproduces 
qualitatively well the previous results obtained both ex- 
perimentally and theoretically. This model seems to 
capture most of the physics of the problem. Because it 
gives a very realistic description of a granular system as 
a three-dimensional packing of hard spheres, it can be a 
quite useful and interesting tool to go beyond a global 
description to a local analysis of the packings' structure 
during the compaction process. 



IV. LOCAL ANALYSIS 

To study locally the packings of spheres more or less 
compacted, we use two kind of descriptions: packing frac- 
tion (or density) profiles are calculated vertically and 
transversly to the box, and size and volume distributions 
of the pores in a packing are evaluated and then analysed. 



A. Density profiles 



C. Ageing 

In these kind of systems in slow evolution to a final 
equilibrium, it is possible to demonstrate ageing effects 
by comparing the system at different ages. This compari- 
son can be made by use of temporal correlation functions 
of global values (p,(Z), ...) between the initial packing 



Using the surface packing fraction calculated by stere- 
ological cuts (as in the evaluation of and ( < i )C )), we 
can have access to vertical (horizontal cuts) and trans- 
verse (vertical cuts) density profiles. 
Some examples of vertical profiles are shown in figure |(]. 
These have been obtained with fixed boundary conditions 
(F.B.C.), but the use of periodic boundaries (P.B.C.) in- 
duces no significant differences. The profiles are char- 
acterised in particular by a negative vertical gradient a 
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and by large peaks near the bottom of the box. These 
peaks reflect a partially ordered packing due to the wall 
and are very close to previous experimental observations 
pB| . The gradient can be roughly estimated in an inter- 
mediate zone (5 < z/R < 22 for F.B.C. and 5 < z/R < 
26 for P.B.C.) after smoothing the profile. This gradient, 
directly linked to e, is qualitatively different from pre- 
vious numerical results where a local densification 
is obtained at the interface. It could be objected that 
this gradient comes directly from the modeling of the 
tap through an uniform dilation. Nevertheless, despite 
the fact that a is difficult to estimate very precisely, it 
does not seem to be monotonic with e 1 / 2 , but has more or 
less the same kind of bell-shaped dependance as the other 
steady-state values. This behavior can not be caused only 
by the dilation. But, in contrast to the other average val- 
ues, it seems that a presents no hysteresis effect which 
denotes a relatively different behavior. In conclusion, the 
origin of the anisotropy of the packing, observed through 
this gradient a, is not well understood. It may come from 
both the uniform dilation of the system and its specific 
redeposition under the simulated particle motion under 
gravity. 



Being inspired by a Fermi level profile |26|, we can pro- 
pose the following average fit for a typical vertical profile: 



*00 



$0 — olz 



I + exp({3(z -z*)) 



where z = z/R (7) 



Figure [t] presents a few transverse profiles in fixed 
boundaries; they are qualitatively close to experimen- 
tal profiles |27|| . Here again, some peaks indicate a local 
organisation in layers due to the walls; this effect has 
approximately a three layers range. The average lateral 
density increase (at a distance less than 7R from the walls 
corresponding roughly to this wall effects range) is noted 
^lateral as is the central increase 8<$> cen t r ai- The last 
one is systematically smaller than the other. Both of 
them are calculated in comparison with the initial pro- 
files and reflect the spatial repartition of the bulk com- 
paction. These profiles with periodic boundaries reveal 
no peak, due to the absence of walls. The central zone 
is a bit larger but keeps the same qualitative shape and 
densifies as well during a tapping sequence. This obser- 
vation of an obvious compaction even in P.B.C. ensures 
that compaction is not, or at least not principally, due to 
wall effects. This was not evident considering the small 
aspect ratio used in the experience of the Chicago group. 
Quantitatively, the absolute value of the packing fraction 
is larger in the periodic conditions but its increase due 
to compaction is a bit smaller. 

As global values, 5$ lateral and 5§ cen t ra i have the same 
dependence on e 1 / 2 (bell shaped curves) than the others. 
It is also possible to study their evolution with the num- 
ber of taps. The results, presented in figure ||, point out, 
once again, the nearly frozen dynamics for small values 
ofe 1 ^. 

Moreover, we can remark that the initial packing in 



F.B.C. presents a great metastability. This one is partic- 
ularly noticeable on the transverse profile (see fig. 0) with 
an "under-population" of the spheres near the walls. This 
explains the significant compaction of 5&i a t e rai caused by 
the first tap (see fig. ||) . This metastability is due to the 
construction of the initial packing: it was built by a grav- 
itational algorithm [16] with periodic conditions, and a 
slight agitation was then induced in the packing to adapt 
it to fixed boundary conditions (F.B.C). This last stage 
was not sufficiently efficient. 



B. Size and volume distributions of the pores 

Another way to analyse a packing of particules is to 
study the interstitial voids. This void-space is more dif- 
ficult to apprehend because, in contrast to a particle, a 
cavity has no geometric limit. We then introduce the 
notion of pore as the "void" between four neighboring 
spheres. Previous studies have already been made on 
this issue, both theoretically and experimentally. Go- 
toh |2^] introduced the pore size distribution Po as the 
probability for a randomly positioned sphere of radius 
r' to intercept no particle center. He proposed a theo- 
retical expression for Po derived from the Percus-Yevick 
approximation which agrees well with previous results on 
random close packings {29 30 1 : 



0<ct<1, P (cr) = l-$cr 

1<<t, P (a) = (1 - ^)exp[--^— 



y± - -if, 

+9/2$cr 2 + 1 - 5/2$)] 



(8) 

<!> .(-(i-2$y 3 

(9) 



where a is the ratio r'/R (R is the radius of the hard 
spheres) and $ is the average volume fraction. Figure 
^ confronts this expression with the distribution calcu- 
lated in one of the packing obtained after a compaction 
sequence of 4000 taps ($ = 60.6%) and with periodic 
boundaries. The distribution calculated in the initial 
packing (<I> = 58.4%) is also represented and is quite 
close to the other. This slight difference means that the 
distribution Po is insufficiently sensitive to small struc- 
tural changes as compaction. 

Hence we find that it is more efficient to work with a 
direct statistical analysis on the size of the pores. To do 
this, we use the Voronoi tesselation of a packing |3~]}] : a 
Voronoi polyhedon around a sphere is the region of space 
in which all the points are closer to this given sphere than 
to any other. Two neighbors correspond to two Voronoi 
polyhedra that share a face. Each vertex is equidistant 
from the center of four neighboring spheres and therefore 
constitutes a pore. More precisely, we define the pore as 
the virtual sphere in contact with these four neighboring 
spheres which interpenetrates none of them. The size of 
the pore is then the radius of this "void- sphere" . The 
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volume of this sphere reflects partially the total void- 
volume situated inside the tetrahedron formed by the 
centers of the four neighboring spheres. In a packing, it 
is possible to calculate the size distribution of the pores 
where £ = r/R with r the radius of a pore and R 
the radius of the hard spheres. The normalisation of the 
distribution gives J p^(^)d^ = Np, the total number of 
pores in the measurement volume. This distribution is 
linked to Gotoh's one. Thus, Pq(cf = r'/R) is more or 
less the sum of the pores of size greater than r'. There- 
fore Pq is a cumulative distribution in comparison with 
p^, which is expected to be more sensitive to the local 
structure. A rather similar analysis by use of the size 
distribution of the pores was previously sketched out by 
Barker and Mehta @. 

If we now use directly the normalized volume v — £ 3 = 
4/3iTfl 3 as new var i a ble, the corresponding statistical den- 
sity is p v (v) = fp£(£) (here the "volume" of a pore is 
reduced to 4/3irr 3 ). The distribution vp v (v) reflects the 
contribution of the pores to the total porosity according 
to their size; this last distribution seems to be the more 
pertinent in problems of free volume and compaction. We 
have noted that, by integrating vp v (v) , a new global value 
is obtained and corresponds to the average normalized 
volume of a spherical pore: v — J vp v (v)dv = (v)/V, 
where (v) is the average volume of a pore, V is the vol- 
ume of the hard spheres, and Np is the number of pores. 
The average pore volume v has the same dynamics and 
the same kind of rever sible-irreversible behavior as de- 
scribed in section [II B . Figure a) shows the distribu- 
tions vp v (v) for a given packing at different stages of its 
compaction (in F.B.C. and with an excitation strength 
e = 2 x 1CP 2 ). The statistics are calculated in a smaller 
box of height 18R located at a distance 2R from the 
walls to avoid some boundary problems. In this mea- 
surement volume, Np varies approximately from 9550 to 
10250 for the different packings analysed. These different 
volume distributions vp v (v) are slightly affected by the 
taps in the small pore domain, whereas the variation of 
the packing fraction clearly appears in the progressive 
reduction of the tail of the distribution in the large pore 
zone. So, there is more or less a persistence of the distri- 
bution for the values of v approximately limited by the 
volume of an octohedral site (£o = 1 ~ 0.414 and 

vo ~ 0.0711). The distribution vp v (v) is bell-shaped 
with an over-population for the largest pores (the tail) 
and with a minimum size of the pores corresponding 
to a tetrahedral site ( £t = \/3/2 — 1 « 0.225 and 
vt ~ 0.0114). With respect to the small pores, this is in 
contradiction with a Poisson distribution proposed in a 
previous theorical model for logarithmic dynamics fl32| . 
But, in figure [l0]b), we note that, in the range of volume 
corresponding to the tail of the distribution, p v (v) is 
compatible with a Poisson law: p v (v) oc e^~ v ^ v °\ where 
vq is directly linked to v, the average normalized volume 
already defined, or to (<f> c ), the average packing fraction. 



These results must be compared with previous work 
on the issue. First, Bernal analyzed the arrangements 
of spheres by characterizing the cavities between the 
spheres. To do this, he studied the different polyhedra 
formed by the sphere centers as corner. He found five 
canonical holes. Table 1 153| presents his results obtained 
on a mechanical model of hard spheres and concerning 
the statistical weight (in number and in volume) and 
the £-value corresponding to each hole. In fact, these 
canonical holes are more or less distorted (otherwise 
P(,{€) would be an addition of Dirac peaks) and the £- 
value corresponding to the regular hole is therefore a 
lower limit. So we can note that the tetrahedron and 
the tetragonal dodecahedron correspond to the smallest 
values of £ for which we have shown that the volume 
distribution vp v (v) is not greatly affected by the tapping 
process. In contrast, the octahedron, the trigonal prism 
and the archimedian antiprism (the two last in very low 
proportions) appear in the large pores range and are 
consequently more sensitive to the compaction state of 
the packing. So, according to this Bernal's classification 
on the structure of the pores, we can suggest that com- 
paction is principally due to rearrangements of the three 
largest canonical pores. 



V. CONCLUSIONS AND PERSPECTIVES 

A simple model of hard spheres under vertical taps 
based solely on geometric constraint is sufficient to quali- 
tatively describe previous experimental and numerical re- 
sults: the same kind of compaction dynamics, hysteresis 
effect on the steady-state values, and ageing behaviors. 
The originality of this model, i.e. a realistic description 
of a granular system as a three-dimensional packing of 
hard spheres, permits a structural analysis of the pack- 
ings. A semi- local study on density profiles suggests the 
existence of a negative vertical gradient in the packings 
but with no clear hysteresis effect. It also confirms a 
compaction in the bulk which can not be caused only by 
wall effects, which are particularly noticeable with fixed 
boundaries (F.B.C). A more local analysis on the void 
space of the packings shows, in a model of spherical pores, 
a volume distribution sensitive to the packing fraction for 
the large pores and nearly stationary for the small ones. 
Compaction could then be principally explained by col- 
lective rearrangements of the largest pores. 
Further to this numerical work, an experimental study of 
the compaction induced by vertical taps is being carried 
out. The packing fraction is deduced from a measure 
of absorption of a horizontal gamma-rays beam. In addi- 
tion to the average volume fraction in the bulk, our set-up 
permits the evaluation of the vertical density profile in 
the packing; these measurements would be crucial in or- 
der to test the results of our numerical model, especially 
concerning the existence of a negative vertical gradient. 
Furthermore, compaction is studied in quite different ex- 
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perimental conditions than in the previous work of the 
Chicago group. In this later one, the set-up is a tube 
of height 80 cm and with an approximate diameter of 2 
cm filled with 1 mm diameter soda lime glass beads. So, 
transverse wall effects are very significant and the ver- 
tical pressure on the packing is saturated for almost all 
of the height of the heap, the overload being completely 
held up by the walls. Inversely, the cylinder used in our 
set-up has a diameter of 10 cm and a height of 15 cm; 
about 80 percent of it is filed with 1 mm diameter glass 
beads. Here, the wall effects become negligible and the 
vertical pressure is definitely not saturated in the pack- 
ing. It will be interesting to see to which extent these 
differences can qualitatively and quantitatively affect the 
compaction under vertical tapping. 
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FIG. 1. A typical displacement during the redeposition stage of the algorithm. 



FIG. 2. Bottom packing fraction 3>t, versus t, the number of taps, in logarithmic (up) and linear (down) time scale for three 
excitation rates (e = 5 x 10~ 3 , 5 x 10~ 2 , and 1.5 x 10 _1 ). The solid lines are the simulation results and the dotted lines are 
the inverse-logarithmic fits. 



FIG. 3. The power spectrum of the fluctuations of (Z) versus frequency (the inverse of the number of taps) (up) and the 
simple standard deviation of (Z), o"(Z)' which is nearly linear with s 1 ' 2 (down). 
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FIG. 4. Steady-state values of ($ c ) obtained after 4000 taps with different values of e 1 / 2 (solid black squares) and hysteresis 
during a sequence of increase (open circles), decrease (open up-triangles) and increase again (open down-triangles) of the 
excitation with an increment Ae 1 / 2 every 4000 taps. 

FIG. 5. Ageing effects on the time-correlation function A(t,tw) for several waiting times tw '■ the different curves (up) and 
a collapse according to the fit in dotted line (down) . 

FIG. 6. Two examples of vertical density profiles with their fit : for the initial packing (dotted line) and for a packing 
obtained after 4000 taps with e = 10 _1 (solid line). 

FIG. 7. Two examples of transverse density profiles : for the initial packing (dotted line) and for a packing obtained after 
4000 taps with e = 10 _1 (solid line). The two vertical dotted lines indicate the frontiers in the calculus of 5$i atera i and 

central • 

FIG. 8. Lateral (up) and central (down) packing fraction increases (in comparison with the initial packing) versus t, the 
number of taps for 4 values of e. 

FIG. 9. Theorical pore size distribution Po (thick solid line) compared with the numerical calculation for a packing with the 
same average volume fraction (open circles). No significant difference with the calculations for the initial packing (thin solid 
line) which has a quite smaller volume fraction. 

FIG. 10. a) Volume distribution of the pores vp v {v) for a packing at different stages of its compaction (F.B.C. and 
e = 2 x 10~ 2 ): t+1 = 1, 100, and 10000. As the packing progressively densities (i.e. as t increases), the tail of the dis- 
tribution corresponding to the largest pores tends to vanish (as symbolized by the arrow). Here, T and O indicate the v- values 
for tetrahedral and octahedral sites in a dense packing (F.C.C. or H.C.). b) Plot of ln(p v (v)) versus v for the same packings. 
The different tails are compatible with a Poisson law. 







Number (%) 


Volume (%) 


tetrahedron 


0.225 


73.0 


48.4 


half-octahedron 


0.414 


20.3 


26.9 


trigonal prism 


0.528 


3.2 


7.8 


tetragonal dodecahedron 


0.353 


3.1 


14.8 


archimedian antiprism 


0.645 


0.4 


2.1 



TABLE I. Characteristics of the Bernal canonical holes. 
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